{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "HMM.ipynb",
      "version": "0.3.2",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_Esy4bIi3E4L",
        "colab_type": "text"
      },
      "source": [
        "# HMM"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jN7qEs6b3UAf",
        "colab_type": "text"
      },
      "source": [
        "Hidden Markov Model, HMM 是可用于标注问题的统计学习模型，属于**生成模型**。由**初始概率分布**，**状态转移概率分布**和**观测概率分布**确定。\n",
        "$\\lambda = (A, B, \\pi)$.  \n",
        "\n",
        "\n",
        "### **两个基本假设：**  \n",
        "1). 齐次马尔可夫性假设，即假设隐藏的马尔可夫链在任意时刻$t$的状态只依赖于其前一时刻的状态，与其他时刻状态及观测无关，也与时刻$t$无关：  \n",
        "$P(i_{t}|i_{t-1},o_{t-1},...,i_{1},o_{1}) = P(i_{t}|i_{t-1}),  t = 1,2,...,T$  \n",
        "\n",
        "2）. 观察独立性假设，即假设任意时刻的观测只依赖与该时刻的马尔可夫链的状态，与其他观测及状态无关：  \n",
        "$P(o_{t}|i_{T},o_{T},i_{T-1},o_{T-1},...,i_{t+1},o_{t+1},i_{t},i_{t-1},o_{t-1},...,i_{1},o_{1} = P(o_{t}|i_{t})$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "W1ADpieg1V3d",
        "colab_type": "text"
      },
      "source": [
        "### **三个基本问题**:   \n",
        "\n",
        "1. **概率计算问题**。给定模型$\\lambda = (A, B, \\pi)$ 和观测序列 $O=(o_{1},o_{2},...,o_{T})$, 计算在模型 $\\lambda$ 下观测序列 $O$ 出现的概率 $P(O|\\lambda)$.  \n",
        "Evaluate $P(O|\\lambda)$.\n",
        "\n",
        "2. **学习问题**。已知观测序列 $O=(o_{1},o_{2},...,o_{T})$, 估计模型 $\\lambda = (A, B, \\pi)$ 的参数，使得在该模型下观测序列概率 $P(O|\\lambda)$ 最大。即用极大似然估计的方法估计参数。  \n",
        "$\\lambda_{MLE} = argmax_{\\lambda}P(O|\\lambda)$.\n",
        "\n",
        "3. **预测问题**。也称为解码(decoding) 问题。已知模型 $\\lambda = (A, B, \\pi)$ 和观测序列 $O=(o_{1},o_{2},...,o_{T})$，求给定观测序列条件概率 $P(I|O)$ 最大的状态序列 $I = (i_{1}, i_{2}, i_{3},...,i_{T})$. 即给定观测序列，求最有可能的对应的状态序列。  \n",
        "$argmax_{I}P(I|O,\\lambda)$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oguS4aSW4Y_Q",
        "colab_type": "text"
      },
      "source": [
        "问题1. 前向（forward）和后向（backward）算法。  \n",
        "问题2. Baum-Welch 算法。  \n",
        "问题3. 近似算法，维特比算法。 "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "L6QU6kS17PgN",
        "colab_type": "text"
      },
      "source": [
        "---------------------------------------------------------------------------------------------------------------------------------"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6aNZwSt17DFG",
        "colab_type": "text"
      },
      "source": [
        "以下来自徐亦达的课程  \n",
        "视频地址：https://www.youtube.com/watch?v=Ji6KbkyNmk8  \n",
        "lecture： https://github.com/roboticcam/machine-learning-notes/blob/master/dynamic_model.pdf"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lyPGx2RqGWrp",
        "colab_type": "text"
      },
      "source": [
        "# 概率计算问题"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "F_XYKDZa9wdQ",
        "colab_type": "text"
      },
      "source": [
        "---------------------------------------------------------------------------------------------------------------------------------"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "596GcIzAqwRs",
        "colab_type": "text"
      },
      "source": [
        "### 例 10.2  \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MuSecU91q2ey",
        "colab_type": "text"
      },
      "source": [
        "考虑盒子和球模型$\\lambda = (A, B, \\pi)$ , 状态集合 $Q = {1,2,3}$, 观测集合 $V = {红,白}$,  \n",
        "\n",
        "\n",
        "$A=\\begin{bmatrix}\n",
        " 0.5&  0.2& 0.3\\\\ \n",
        " 0.3&  0.5& 0.2\\\\ \n",
        " 0.2&  0.3& 0.5\n",
        "\\end{bmatrix}, \n",
        "B=\\begin{bmatrix}\n",
        " 0.5&  0.5\\\\ \n",
        " 0.4&  0.6\\\\ \n",
        " 0.7&  0.3\n",
        "\\end{bmatrix},\n",
        "\\pi=\\begin{bmatrix}\n",
        " 0.2\\\\ \n",
        " 0.4\\\\ \n",
        " 0.4\n",
        "\\end{bmatrix}$  \n",
        "\n",
        "设$T=3, Q={红,白,红}$,试用前向算法计算$P(O|\\lambda)$.  \n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gDkjQ0uJy_-M",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "A = [[0.5, 0.2, 0.3], [0.3, 0.5, 0.2],[0.2, 0.3, 0.5]]\n",
        "B = [[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]]\n",
        "pi = [0.2, 0.4, 0.4]\n",
        "Q = [0,1,0]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zE5Yt4jdu_VP",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZwRR7eFrHCnJ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "class HMM_fw:\n",
        "    def __init__(self, A, B, pi):\n",
        "        self.A = A    # 状态转移概率\n",
        "        self.B = B    # 观测概率\n",
        "        self.pi = pi  # 初始状态\n",
        "        \n",
        "    def forward(self, Q):\n",
        "        T = len(Q)   # 观测序列长度，时刻T\n",
        "        M = len(self.A)   # 状态数\n",
        "        alpha = np.zeros((T, M))\n",
        "        \n",
        "        for t in range(T):\n",
        "            for m in range(M):\n",
        "                if t == 0:\n",
        "                    alpha[t][m] = self.pi[m] * self.B[m][Q[t]]\n",
        "                    print(\"alpha[{}][{}] = pi[{}] * B[{}](Q{}) = {:.2f}\".format(t+1, m+1, m+1, m+1, Q[t]+1, alpha[t][m]))\n",
        "                else:\n",
        "                    alpha[t][m] = sum([alpha[t-1][i] * self.A[i][m] for i in range(len(alpha[t-1]))]) * self.B[m][Q[t]]\n",
        "                    print(\"alpha[{}][{}] = {:.5f}\".format(t+1, m+1, alpha[t][m]))\n",
        "                    \n",
        "        p = sum(alpha[T-1])\n",
        "        #print(p)\n",
        "        return p"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HfqXwAO1y9G2",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 197
        },
        "outputId": "a6dcb119-4da5-4dd6-a18c-593bfd70eaa4"
      },
      "source": [
        "m = HMM_fw(A, B, pi)\n",
        "m.forward(Q)"
      ],
      "execution_count": 79,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "alpha[1][1] = pi[1] * B[1](Q1) = 0.10\n",
            "alpha[1][2] = pi[2] * B[2](Q1) = 0.16\n",
            "alpha[1][3] = pi[3] * B[3](Q1) = 0.28\n",
            "alpha[2][1] = 0.07700\n",
            "alpha[2][2] = 0.11040\n",
            "alpha[2][3] = 0.06060\n",
            "alpha[3][1] = 0.04187\n",
            "alpha[3][2] = 0.03551\n",
            "alpha[3][3] = 0.05284\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0.130218"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 79
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9LcfAGUVzaU3",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "class HMM_bw:\n",
        "    def __init__(self, A, B, pi):\n",
        "        self.A = A\n",
        "        self.B = B\n",
        "        self.pi = pi\n",
        "        \n",
        "    def backward(self, Q):\n",
        "        T = len(Q)   # 观测序列长度，时刻T\n",
        "        N = len(self.A)   # 状态数\n",
        "        beta = np.zeros((T, N))\n",
        "        \n",
        "        for t in range(T-1, -1, -1):\n",
        "            for n in range(N):\n",
        "                if t == T - 1:\n",
        "                    beta[t][n] = 1\n",
        "                    print(\"beta[{}][{}] = {:.2f}\".format(t+1, n, beta[t][n]))\n",
        "                else:\n",
        "                    beta[t][n] = sum(self.A[n][j] * self.B[j][Q[t+1]] * beta[t+1][j] for j in range(N))\n",
        "                    print(\"beta[{}][{}] = {:.5f}\".format(t+1, n, beta[t][n]))\n",
        "                    \n",
        "        p = sum(self.pi[i] * self.B[i][Q[0]] * beta[0][i] for i in range(N))\n",
        "        #print(p)\n",
        "        return p"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rIyCzvnbBD8K",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 197
        },
        "outputId": "48682810-a2c9-47ca-d46c-2cb3eeb545b5"
      },
      "source": [
        "m = HMM_bw(A, B, pi)\n",
        "m.backward(Q)"
      ],
      "execution_count": 77,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "beta[3][0] = 1.00\n",
            "beta[3][1] = 1.00\n",
            "beta[3][2] = 1.00\n",
            "beta[2][0] = 0.54000\n",
            "beta[2][1] = 0.49000\n",
            "beta[2][2] = 0.57000\n",
            "beta[1][0] = 0.24510\n",
            "beta[1][1] = 0.26220\n",
            "beta[1][2] = 0.22770\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0.130218"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 77
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uc5CfmwRGdb6",
        "colab_type": "text"
      },
      "source": [
        "# 学习问题"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qTYQ6fO8GfMI",
        "colab_type": "text"
      },
      "source": [
        "HMM的学习，根据训练数据是否包括观测序列和对应的状态序列还是只有观测序列，可以分别为**有监督学习**和**无监督学习**来实现。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qTRkuN4GIfnW",
        "colab_type": "text"
      },
      "source": [
        "### 监督学习  \n",
        "\n",
        "假设已给训练数据包含$S$个长度相同的观测序列和对应的状态序列${(O_{1}, I_{1}), (O_{2}, I_{2}),..., (O_{S}, I_{S})}$. 那么可以利用**极大似然估计**法来估计HMM的参数。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "neoLTatUJgks",
        "colab_type": "text"
      },
      "source": [
        "### 无监督学习  \n",
        "\n",
        "假设已给训练数据只包含$S$个长度为$T$的观测序列${O_{1}, O_{2},..., O_{S}}$, 而没有对应的状态序列， 目标是学习HMM $\\lambda = (A, B, \\pi)$ 的参数。 我们将观测序列数据看作观测数据$Q$, 状态序列数据看作不可观测的隐数据$I$, 那么HMM则是一个含有隐变量的概率模型：  \n",
        "\n",
        "$P(O|\\lambda) = \\sum_{I}P(O|I, \\lambda)P(I|\\lambda)$  \n",
        "\n",
        "他的参数可以由EM算法来学习。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9862N0kFBFBb",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0lKwVVTPPpjt",
        "colab_type": "text"
      },
      "source": [
        "# 预测问题\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "796IZ5QqPrdS",
        "colab_type": "text"
      },
      "source": [
        "### 近似算法  \n",
        "近似算法的想法是， 在每个时刻$t$ 选择在该时刻最有可能出现的状态$i^{*}_{t}$，从而得到一个状态序列 $I^{*} = (i^{*}_{1}, i^{*}_{2}, ..., i^{*}_{T})$, 将他作为预测的结果。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G5XHAHOdQlKE",
        "colab_type": "text"
      },
      "source": [
        "### 维特比算法  \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TlCSEMGhQnRM",
        "colab_type": "text"
      },
      "source": [
        "维特比算法实际是用动态规划，解HMM的预测问题，即用动态规划求概率最大路径。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gSOX42cFSM2l",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "class HMM_viterbi:\n",
        "    def __init__(self, A, B, pi):\n",
        "        self.A = A\n",
        "        self.B = B\n",
        "        self.pi = pi\n",
        "        \n",
        "    def viterbi(self, Q):\n",
        "        T = len(Q)   # 观测序列长度，时刻T\n",
        "        N = len(self.A)   # 状态数\n",
        "        sigma = np.zeros((T, N))\n",
        "        delta = np.zeros((T, N))\n",
        "        for t in range(T):\n",
        "            for n in range(N):\n",
        "                if t == 0:\n",
        "                    sigma[t][n] = self.pi[n] * self.B[n][Q[t]]\n",
        "                    delta[t][n] = 0\n",
        "                    print(\"sigmia[{}][{}] = {:.2f}\".format(t+1, n+1, sigma[t][n]))\n",
        "                    print(\"delta[{}][{}] = {}\".format(t+1, n+1, delta[t][n]))\n",
        "                    \n",
        "                else:\n",
        "                    sigma[t][n] = np.max([sigma[t-1][j] * self.A[j][n] for j in range(N)]) * self.B[n][Q[t]]\n",
        "                    print(\"sigma[{}][{}] = {:.5f}\".format(t+1, n+1, sigma[t][n]))\n",
        "                    \n",
        "                    delta[t][n] = np.argmax([sigma[t-1][j] * self.A[j][n] for j in range(N)]) + 1\n",
        "                    print(\"delta[{}][{}] = {}\".format(t+1, n+1, delta[t][n]))\n",
        "                    \n",
        "        P = np.max(sigma[T-1])\n",
        "        print(P)\n",
        "        pth = []\n",
        "        for t in range(T-1, -1, -1):\n",
        "            if t == T - 1:\n",
        "                i_t = np.argmax(sigma[t])\n",
        "                pth.append(i_t + 1)\n",
        "            else:\n",
        "                i_t = int(delta[t+1][i_t]) - 1\n",
        "                pth.append(i_t + 1)\n",
        "        \n",
        "        return pth"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LDEr3xz9SNSE",
        "colab_type": "text"
      },
      "source": [
        "#### 例 10.3"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oFJgnsvySRNF",
        "colab_type": "text"
      },
      "source": [
        "$A=\\begin{bmatrix}\n",
        " 0.5&  0.2& 0.3\\\\ \n",
        " 0.3&  0.5& 0.2\\\\ \n",
        " 0.2&  0.3& 0.5\n",
        "\\end{bmatrix}, \n",
        "B=\\begin{bmatrix}\n",
        " 0.5&  0.5\\\\ \n",
        " 0.4&  0.6\\\\ \n",
        " 0.7&  0.3\n",
        "\\end{bmatrix},\n",
        "\\pi=\\begin{bmatrix}\n",
        " 0.2\\\\ \n",
        " 0.4\\\\ \n",
        " 0.4\n",
        "\\end{bmatrix}$  \n",
        "\n",
        "已知观测序列$O=(红, 白, 红)$,试求最优状态序列，即最优路径 $I^{*}=(i^{*}_{1}, i^{*}_{2}, i^{*}_{3})$."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MXWapkQESxlq",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "A = [[0.5, 0.2, 0.3], [0.3, 0.5, 0.2],[0.2, 0.3, 0.5]]\n",
        "B = [[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]]\n",
        "pi = [0.2, 0.4, 0.4]\n",
        "Q = [0,1,0]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "a55NPE7AS-Eu",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 377
        },
        "outputId": "5f86553f-345a-466d-fc8a-755b65a1d167"
      },
      "source": [
        "m = HMM_viterbi(A, B, pi)\n",
        "m.viterbi(Q)"
      ],
      "execution_count": 135,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "sigmia[1][1] = 0.10\n",
            "delta[1][1] = 0.0\n",
            "sigmia[1][2] = 0.16\n",
            "delta[1][2] = 0.0\n",
            "sigmia[1][3] = 0.28\n",
            "delta[1][3] = 0.0\n",
            "sigma[2][1] = 0.02800\n",
            "delta[2][1] = 3.0\n",
            "sigma[2][2] = 0.05040\n",
            "delta[2][2] = 3.0\n",
            "sigma[2][3] = 0.04200\n",
            "delta[2][3] = 3.0\n",
            "sigma[3][1] = 0.00756\n",
            "delta[3][1] = 2.0\n",
            "sigma[3][2] = 0.01008\n",
            "delta[3][2] = 2.0\n",
            "sigma[3][3] = 0.01470\n",
            "delta[3][3] = 3.0\n",
            "0.014699999999999998\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[3, 3, 3]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 135
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GKPLvAFzXuva",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}